% 实线性Klein-Gordon方程的1D显式求解。
% 在一定条件下显式求解应该是稳定的，参考 doi: 10.1007/s00211-011-0411-2
% 当粒子初始的位置“足够准确时”(\Delta x足够小)，有可能观察到波前超光速运行（根据不确定性原理，此时\Delta p 足够大，即粒子可能具有非常高的动量），而这是违反相对论的
% 因此一般认为“单粒子”模型在相对论性量子力学中并不是合理的

% Gitee Repo

clc
clear

L=10;
dx = 0.01;
dt = 0.001;

m=1;
hbar = 1;
c=1;

sigma = (c*dt/dx)^2;
dtmch = dt^2*m^2*c^4/hbar^2;

x = (-L:dx:L)';
n = size(x,1);

u0 = zeros(n,1);
u1 = zeros(n,1);
u2 = zeros(n,1);

u0 = e.^(-100*x.^2);
u0(abs(u0)<0.001) = 0;
u1 = u0;
u2 = u0;

idxr = find(u0 ~= 0, 1, 'last');
idxl = find(u0 ~= 0, 1, 'first');
wavefrontFound=1;
if isempty(idxr) || isempty(idxl)
    warning('wavefront not found')
    wavefrontFound=0;
end

figure()
imgindex = 0;

TICK=10000;
for tick = 3:TICK

    u2(2:n-1) = 2*u1(2:n-1) - u0(2:n-1) + sigma*(u1(1:n-2)-2*u1(2:n-1) + u1(3:n)) - dtmch*u1(2:n-1);

    u2(1) = 0;
    u2(n) = 0;

    u0 = u1;
    u1 = u2;

    if tick==3 || mod(tick,100) == 0
        clf
        hold on
        axis equal
        axis([-L L -2 2])

        plot(x,u2)

        line([-L L],[0 0])
        if wavefrontFound == 1
            scatter(x(idxl)-c*tick*dt,0);
            scatter(x(idxr)+c*tick*dt,0);
        end

        drawnow
        pause(0.1)
        imgindex++;
    end
%    break
end

